;; Tests of the non-robust linear algebra functions
;; Copyright (C) 2025  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.linalg.sole-test
  (:require
   [clojure.test :refer [deftest is testing]]
   [mazama.linalg.sole :as sole]
   [mazama.linalg.utilities :as utilities]))

(def EPSILON 1/1000000)

(defn nearly==
  ""
  ([a b] (nearly== a b EPSILON))
  ([a b epsilon]
   (< (- epsilon) (- a b) epsilon)))

(defn- nullspace-vector?
  "Is the vector `v` in the nullspace of the matrix `A`?"
  [A v]
  (every? zero? (utilities/*mv A v)))

(defn- check-Ax==b [A x b]
  (let [Ax (utilities/*mm A x)]
    (is (== (:ncols Ax) (:ncols b)))
    (doseq [[x y] (map vector (:data Ax) (:data b))]
      (is (== x y)))))

(defn- check-general-solution
  [A b expected-nullspace-dimension]
  (let [x (sole/solve A b)]
    (check-Ax==b A x b)
    (is (== 3 (count (keys x))))
    (is (:ncols x))
    (is (:data x))
    (is (:nullspace x))
    (let [nullspace (:nullspace x)]
      (is (== expected-nullspace-dimension (count nullspace)))
      (is (every? #(nullspace-vector? A %) nullspace))
      (is (utilities/linearly-independent? nullspace)))))

(defn- check-unique-solution
  [A b]
  (let [x (sole/solve A b)]
    (check-Ax==b A x b)
    (is (== 2 (count (keys x))))
    (is (:ncols x))
    (is (:data x))))

(deftest sole-no-solutions
  (testing "lem.ma The Simplest Possible System (without a solution)"
    (let [A {:ncols 1 :data [0]}
          b {:ncols 1 :data [10]}]
      (is (nil? (sole/solve A b)))))

  (testing "http://sites.science.oregonstate.edu/math/home/programs/undergrad/CalculusQuestStudyGuides/vcalc/gauss/gauss.html"
    (let [A {:ncols 3 :data [ 2 -1  1
                              3  2 -4
                             -6  3 -3]}
          b {:ncols 1 :data [1 4 2]}]
      (is (nil? (sole/solve A b)))))

  (testing "https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/linear-systems/gaussian-elimination"
    (let [A {:ncols 3 :data [ 1  1 -3
                              2  1 -1
                              3  2 -4]}
          b {:ncols 1 :data [4 2 7]}]
      (is (nil? (sole/solve A b))))))

(deftest sole-one-solution
  (testing "More equations than unknowns"
    (let [A {:ncols 2 :data  [4  7
                              5  8
                              6  9]}
          b {:ncols 1 :data [1 2 3]}]
      (check-unique-solution A b)))

  (testing "https://ocw.mit.edu/ans7870/18/18.013a/textbook/HTML/chapter32/section01.html"
    (let [A {:ncols 3 :data [ 2  3  1
                              3 -1 -2
                             -1 -1  1]}
          b {:ncols 1 :data [7 4 1]}]
      (check-unique-solution A b)))

  (testing "https://ocw.mit.edu/ans7870/18/18.013a/textbook/HTML/chapter32/exercise02.html"
    (let [A {:ncols 3 :data [ 1  3 -1
                              3  1 -2
                             -1 -1  1]}
          b {:ncols 1 :data [7 4 1]}]
      (check-unique-solution A b)))
  
  (testing "https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/linear-systems/gaussian-elimination"
    (let [A {:ncols 2 :data [1  1
                             3 -2]}
          b {:ncols 1 :data [3 4]}]
      (check-unique-solution A b)))

  (testing "https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/linear-systems/gaussian-elimination"
    (let [A {:ncols 3 :data (list 1 -2  1
                                  2  1 -3
                                  4 -7  1)}
          b {:ncols 1 :data (list 0 5 -1)}]
      (check-unique-solution A b)))

  (testing "https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/linear-systems/gaussian-elimination"
    (let [A {:ncols 3 :data (seq [ 2 -2  0
                                   1 -1  1
                                   0  3 -2])}
          b {:ncols 1 :data [-6 1 -5]}]
      (check-unique-solution A b)))

  (testing "Example 5 from https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/linear-systems/gaussian-elimination"
    (let [A {:ncols 3 :data [ 1/4  1/2  1
                                1    1  1
                                4    2  1]}
          b {:ncols 1 :data (seq [23/4 7 2])}]
      (check-unique-solution A b)))

  (testing "Some extraneous zero rows, using the previous test data"
    (let [A {:ncols 3 :data [   0    0  0
                              1/4  1/2  1
                                0    0  0
                                0    0  0
                                1    1  1
                                0    0  0
                                4    2  1
                                0    0  0
                                0    0  0]}
          b {:ncols 1 :data [0 23/4 0 0 7 0 2 0 0]}]
      (check-unique-solution A b)))

  (let [seed 1
        prng (utilities/prng-min-max seed -10 100)]
    (testing (str "random:vector:" seed)
      (let [[A _ b] (utilities/random-sole-unique-solution prng 13 1)]
        (check-unique-solution A b))))

  (let [seed 2
        prng (utilities/prng-min-max seed -5 5)]
    (testing (str "random:vector:" seed)
      (let [[A _ b] (utilities/random-sole-unique-solution prng 49 1)]
        (check-unique-solution A b))))

  (let [seed 3
        prng (utilities/prng-min-max seed -99999 99999)]
    (testing (str "random:matrix:" seed)
      (let [[A _ b] (utilities/random-sole-unique-solution prng 49 2)]
        (check-unique-solution A b)))))

(deftest sole-invert-matrix

  (testing "https://ocw.mit.edu/ans7870/18/18.013a/textbook/HTML/chapter32/section03.html"
    (let [A {:ncols 2 :data [ 1  3
                              5  6]}
          b {:ncols 2 :data [ 1  0
                              0  1]}]
      (check-unique-solution A b)))

  (testing "https://people.richland.edu/james/lecture/m116/matrices/inverses.html"
    (let [A {:ncols 3 :data [ 3  2 -5
                              1 -3  2
                              5 -1  4]}
          b {:ncols 3 :data [ 1  0  0
                              0  1  0
                              0  0  1]}]
      (check-unique-solution A b))))

(deftest sole-infinite-solutions
  (testing "lem.ma The Simplest Possible System (with infinite solutions)"
    (let [A {:ncols 1 :data [0]}
          b {:ncols 1 :data [0]}]
      (check-general-solution A b 1)))

  (testing "all zeros"
    (let [A {:ncols 2 :data [0 0
                             0 0]}
          b {:ncols 1 :data [0 0]}]
      (check-general-solution A b 2)))

  (testing "some zero columns"
    (let [A {:ncols 4 :data [1 0 2 0
                             2 0 4 0
                             3 0 3 0]}
          b {:ncols 1 :data [3 6 6]}]
      (check-general-solution A b 2)))

  (testing "Column with real and integer zeros"
    (let [A {:ncols 3 :data [  0  2  3
                             0.0  4 -1
                               0  3  0]}
          b {:ncols 1 :data [13 5 6]}]
      (check-general-solution A b 1)))

  (testing "lem.ma Our First Linear System: A Thorough Discussion"
    (let [A {:ncols 4 :data [1  0  4  7
                             0  1  5 11]}
          b {:ncols 1 :data [8 11]}]
      (check-general-solution A b 2)))

  (testing "lem.ma A Wide Linear System"
    (let [A {:ncols 5 :data [2 5 6 7 241]}
          b {:ncols 1 :data [6]}]
      (check-general-solution A b 4)))

  (testing "lem.ma A Look Ahead - An Illustration of Matrix Algebra"
    (let [A {:ncols 5 :data [1 2 3 4 5 3 7 10 13 16]}
          b {:ncols 1 :data [-4 -16]}]
      (check-general-solution A b 3)))

  (testing "https://math.colorado.edu/~nita/SystemsofLinearEquations.pdf Example 0.17 (look at the equations, not the matrices, which have errors."
    (let [A {:ncols 5 :data [2 3 1 4 -9
                             1 1 1 1 -3
                             1 1 1 2 -5
                             2 2 2 3 -8]}
          b {:ncols 1 :data [17 6 8 14]}]
      (check-general-solution A b 2)))

  (testing "https://www.math.ust.hk/~mabfchen/Math111/Week1-3.pdf Example 2.4"
    (let [A {:ncols 6 :data [ 0  0  1 -1  2  1
                              3  6  0  3 -3  2
                              1  2  0  1 -1  0
                              2  4 -2  4 -6 -5]}
          b {:ncols 1 :data [0 7 1 -4]}]
      (check-general-solution A b 3)))

  (testing "https://www.math.ucdavis.edu/~linear/linear-guest.pdf, section 2.1.4, Example 20"
    (let [A {:ncols 4 :data [ 1  0  7  0
                              0  1  3  4
                              0  0  0  0
                              0  0  0  0]}
          b {:ncols 1 :data [4 1 0 0]}]
      (check-general-solution A b 2)))

  (testing "https://math.vanderbilt.edu/sapirmv/msapir/gauss.html"
    (let [A {:ncols 7 :data [ 4  8 -2 -5  0 -23  2
                              1  2  2  5 -5  13 -2
                              2  4 -2 -5  2 -19  2
                              3  6  1  0  0  -1  0
                              4  8  1  2 -5  -1 -1]}
          b {:ncols 1 :data [37 3 21 22 30]}]
      (check-general-solution A b 3)))

  (testing "https://www.math.ust.hk/~mabfchen/Math111/Week1-3.pdf Example 2.2"
    (let [A {:ncols 4 :data [ 1 -1  1 -1
                              1 -1  1  1
                              4 -4  4  0
                             -2  2 -2  1]}
          b {:ncols 1 :data [2 0 4 -3]}]
      (check-general-solution A b 2))))

(deftest sole-rounding-errors

  (testing "https://www.dam.brown.edu/people/alcyew/handouts/linsys2.pdf"
    (let [A {:ncols 3 :data [10  -7       0
                             -3  2.09999  6
                              5  -1       5]}
          b {:ncols 1 :data [7  3.90001 6]}]
      (doseq [[x y] (map vector [0 -1 1] (:data (sole/solve A b)))]
        (is (nearly== x y)))))

  (testing "https://www.dam.brown.edu/people/alcyew/handouts/linsys2.pdf (as rationals)"
    (let [A {:ncols 3 :data [10  -7            0
                             -3  209999/100000 6
                              5  -1            5]}
          b {:ncols 1 :data [7  390001/100000 6]}]
      (is (= [0 -1 1] (:data (sole/solve A b))))))

  (testing "http://site.iugaza.edu.ps/emasry/files/2013/01/Ch9-Gauss_Elimination4.pdf"
    (let [A {:ncols 2 :data [ 2.5  6.2
                              4.8 -8.6]}
          b {:ncols 1 :data [3 5.5]}]
      (doseq [[x y] (map vector [1.168552478 0.01268045242] (:data (sole/solve A b)))]
        (is (nearly== x y)))))

  (testing "http://site.iugaza.edu.ps/emasry/files/2013/01/Ch9-Gauss_Elimination4.pdf"
    (let [A {:ncols 3 :data [ 3.0 -0.1 -0.2
                              0.1  7.0 -0.3
                              0.3 -0.2 10.0]}
          b {:ncols 1 :data [7.85 -19.3 71.4]}]
      (doseq [[x y] (map vector [3.0 -2.5 7.0] (:data (sole/solve A b)))]
        (is (nearly== x y)))))

  (testing "http://site.iugaza.edu.ps/emasry/files/2013/01/Ch9-Gauss_Elimination4.pdf (as rationals)"
    (let [A {:ncols 3 :data [ 30/10 -01/10 -02/10
                              01/10  70/10 -03/10
                              03/10 -02/10 100/10]}
          b {:ncols 1 :data [785/100 -193/10 714/10]}]
      (is (= [30/10 -25/10 70/10] (:data (sole/solve A b)))))))
